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Abstract. Ab mitio calculations and a direct method have been applied to derive the phonon dispersion 
curves and phonon density of states for the TiC crystal. The results are compared and found to be in a good 
agreement with the experimental neutron scattering data. The force constants have been determined from 
the Hellmann-Feynman forces induced by atomic displacements in the 2x2x2 supercell. The calculated 
phonon density of states suggests that vibrations of Ti atoms form acoustic branches, whereas the motion 
of C atoms is confined to optic branches. The elastic constants have been found using the deformation 
method and compared with the results obtained from acoustic phonon slopes. 

PACS. 63.20.-e Phonons and vibrations in crystal lattices - 71.15.N Total energy calculations 
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The transition-metal carbide compounds, of which TiC 
is a representative, are of big scientific and technologi- 
cal interest because of their striking mechanical proper- 
ties, extreme hardness combined additionally with metal- 
lic electrical and thermal conductivity. They have usually 
a rock-salt crystal structure - like some ionic crystals - 
while having properties typical for materials with covalent 
bonding. It suggests that the bonding have some mixed 
nature, which makes these materials very interesting as 
an object of study. 

Over a large range of fractional carbon content from 
TiC to TiCo.5 [B the TiC is stable in the NaCl structure, 
with space group Fm3m. In practice, it is usually non- 
stoichiometric and contains carbon vacancy defects. The 
phonon dispersion curves of nearly stoichiometric TiCo.95 
have been measured along main symmetry directions by 
Pintschovius et al [5] using the inelastic neutron scattering 
technique. In the same paper, the influence of the carbon 
content on the phonon dispersion curves has been addi- 
tionally checked by measuring a sample with lower carbon 
concentration TiCo.89- The results shown that the optic 
modes could differ by 3.5%, while the acoustic modes are 
less sensitive to carbon content. The elastic constants de- 
rived from the measured acoustic dispersion curves agree 
very well with those obtained by ultrasonic measurements 
for TiCo.91 crystal [3]. 

Very interesting properties and simple structure of TiC 
caused that efforts to understand its properties using first- 
principle calculations have been undertaken several times 
[HEllllIZ]. The electronic structure, bulk modulus, and elas- 
tic constants have been found [8] by means of accurate 
first-principles total-energy calculations using the full po- 
tential linear muffin-tin orbital method. The calculated 
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values are generally in very good agreement with experi- 
ment. 

In this paper we intend to extend first-principle cal- 
culations to describe the phonon dispersion curves and 
phonon density of TiC. The method which was is based 
on the total energy calculation and Hellmann-Feynman 
(HF) forces. Dispersion relations are calculated by direct 
method pi fTUllll|[T^ll3) . in which the force constants of 
the dynamical matrix are calculated from HF forces. Elas- 
tic constants and bulk modulus have been estimated by 
straight evaluation of energy derivatives with respect to 
deformation. 

The energies and forces of the TiC crystal were cal- 
culated by the method of total energy minimization, us- 
ing a norm-conserving pseudopotentials as an approxima- 
tion of the atomic core - valence electron interaction [TH 
[T5] . For an excellent review of ab initio total-energy cal- 
culations see Ref.[16 and further references given there. 
For the lattice dynamics calculations the 2x2x2 supercell 
with periodic boundary conditions and 64 atoms has been 
used. For optimizing the structure and for direct calcu- 
lation of energy deformations 1x1x1 supercell has been 
utilized. The ab initio total energy calculations were done 
with the CASTEP package jlj and standard pseudopo- 
tentials provided within this package. These pseudopoten- 
tials were constructed within LDA approximation. The Ti 
pseudopotential treats 3p electrons as belonging to the 
core. We have employed non-local variant of those poten- 
tials parametrized in the reciprocal space with 900 eV and 
600 eV cut-off energy for 1x1x1 and 2x2x2 supercell, 
respectively. Tests which had been made with the cut-off 
energy of 600 eV applied to the 1x1x1 supercell, showed 
that in this case the cohesive energy is only 0.2 eV higher, 
and the lattice constant is longer by 0.0006A(0.01%), com- 
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paring to the 900 eV cut-off. Therefore, we have used 600 
eV cut-off energy for a larger supercell calculation. 

A generalized gradient approximation (GGA) has been 
used for the exchange energy term of the valence states 
of the Hamiltonian [r5|. Although this procedure is not 
quite consistent, there are indications [19] that LDA gen- 
erated pseudopotentials used with GGA exchange term 
in crystals are sufficiently accurate for energies and struc- 
tures. The integration over the Brillouin zone has been 
performed with weighted summation over wave vectors 
generated by Monkhorst-Pack scheme i2Qj using 0.1 A-i 
and 0.06 grid sizes which leads to 4 and 32 wave 
vectors for 1x1x1 supercell and to 1 and 4 wave vectors 
for 2x2x2 supercell, respectively. No significant difference 
due to different number of wave vectors, 1 or 4, has been 
found in calculated phonon frequencies. The results of dis- 
persion curve presented here are generated from the 0.06 
data set with 4 wave vectors. 

Since TiC is a weakly metallic compound it could have 
been appropriate to use a smearing of electron levels in 
the Brillouin zone integration. This procedure is compu- 
tationally more expensive. Thus we have performed a se- 
ries of tests to compare lattice constants, bulk moduli and 
phonons in F and X points obtained using gaussian smear- 
ing in the range from 4 to 0.1 eV and denser k-space grids 
(4, 14, 32, 63 k-points) for one unit cell configuration. The 
comparison of results shown in the Table [1] indicate that, 
for sufficiently dense k-space grid used, these quantities 
are not very sensitive (specially the lattice constant and 
bulk modulus) to the choice of method of Brillouin zone 
integration. Thus we have decided to use much simpler 
and faster non- metallic approach in the final 2x2x2 su- 
percell calculation. Note also that results obtained in the 
metallic regime lead to larger divergency of phonon fre- 
quencies from experimental data. 

The minimization of the total energy with respect to 
the lattice constant was performed within the CASTEP 
minimizer module with 1x1x1 supercell and cross-checked 
with 2x2x2 supercell. The equilibrium lattice constant is 
a = 4.3448A which could be compared with the exper- 
imental value a — 4.3269A[2T]. Such a good agreament 
may be coincidental, since inclusion of 3p electrons in the 
valence band as well as gradient corrections in core states 
tends to change the equilibrium lattice constant [2^1^ 
[24]. 

The HF forces are defined as 

Y,{n,v)^~dEtotld-R,{n,v) (1) 

where n, v are the indices of the unit cell and atom within 
the unit cell, respectively, and i is the Cartesian compo- 
nent. At extremum all HF forces vanish. Non-zeroth HF 
forces arise, when a single atom (m, ji) is displaced by 
itj(m,/i) from its equilibrium position. The arising forces 
are related with the cummulant force constants <?y- by a 
relation [TT I fT ^ 

F, (n, I/) = - ^ <P,.j{xi.,v;xa., n)uj (m, n) (2) 



Cummulant force constants appear as a result of peri- 
odic boundary conditions imposed on the supercell. To 
calculate HF forces an atomic configuration with a sin- 
gle displaced atom must be minimized with respect to the 
electronic part only. Each of such runs provides 3n HF 
forces, where n = 8 or 64 is the number of atoms in the 
1x1x1 or 2x2x2 supercell, respectively. Two runs with 
displacements along z direction, one for Ti and second for 
C were performed. To minimize the systematic errors two 
other runs with negative displacements were carried on. 
The displacement amplitudes were set to 1% of the lattice 
constant. This amplitude has been selected after number 
of careful tests in which displacements ware varied from 
0.1% to 1.6%. It has been proven that the anharmonic 
contributions are negligable up to 1% of displacement. At 
smaller displacements, close to 0.1% the HF forces became 
too small and the uncertainty of force constant evaluation 
is not acceptable. 
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Fig. 1. Absolute eigenvalues of the force constant matrices of 
the TiC crystal as a function of interatomic distance. 

The data of 4 x 3n HF forces Fi(n, z^) and four dis- 
placements Mj(m, fi) form an overdetermined set of equa- 
tions, Eq.®, for the force constants. This system was 
solved by the singular value decomposition algorithm [TTl 
125] ■ which automatically provides the least-squares solu- 
tion. For 2x2x2 supercell the 768 components of the HF 
forces lead to 33 non-zero independent parameters of the 
force constants. Knowledge of these force constants al- 
lows one to test the range of interaction potential. For 
that, similarly to procedure used in [T^, we have calcu- 
lated the absolute values of the eigenvalues of the 3x3 

(n, z/; m, /x) matrices and plotted them against a dis- 
tance between atoms (n, z/) and (m, /i). The result for 
2x2x2 supercell is depicted in Fig.[T] It shows that the in- 
teraction drops down at least by two orders of magnitudes 
within the size of the supercell. The scatter of points is too 
large to assign any rule which could govern the decaying of 
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Quantity 


Exper. 


N on- metallic 


Metallic 


#k 




4 


4 


14 


32 


63 


ao 


4.327 


4.345 


4.35 


4.342 


4.345 


4.342 


B 


2.4 


2.42 


2.41 


3.527 


2.418 


2.427 


ujr 


16.3 


16.5 


19.0 


11.8 


15.6 


15.6 


'^x 


8.5 


8.5 


10.2 


5.1 


9.2 


8.5 


2 

UJx 


10.8 


10.0 


11.1 


10.3 


10.7 


10.7 




16.5 


17.0 


18.2 


12.8 


15.7 


15.3 


4 


18.5 


17.8 


19.3 


16.7 


18.1 


17.7 



Table 1. Comparison of results from metallic and non-metallic approach 



the force constants parameters with distance. From Fig. [T] 
it is clear that the 1x1x1 supercell is too small to be used 
for calculation of dispersion curves since the considered 
force constants terminate at radius of a\/3/2 = 0.866a. 

The knowledge of the cummulant force constants al- 
lows one to define an approximate dynamical matrix, in 
which the summation over all atoms is confined to those 
atoms which reside within the volume of the supercell. 
The approximate dynamical matrix becomes equal to the 
conventional one at discrete wave vectors k/, given be the 
equation exp(27rikL • L) = 1 [TT], where L — (Li,L2,L3) 
are the lattice vectors of the supercell. At the wave vectors 
ki, the phonon frequencies w^(kL), calculated by diago- 
nalization of the approximate dynamical matrix, are the 
same as those calculated from the exact dynamical ma- 
trix. For 2x2x2 supercell the exact wave vectors are at 
r, X, L, W, and two other wave vector points, namely, 
the midpoint between F and X along (1, 0, 0) and (1, 1, 1) 
directions. The advantage of the above described direct 
method is that it does not impose any limit to the range 
of interaction. When the supercell size is smaller than the 
range of interaction, the direct method interpolates the 
dispersion curves between the exact points. 
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Fig. 2. Phonon dispersion relations of TiC crystal calculated 
from 2x2x2 supercell (full line). The same dispersion curves 
divided by a factor 1.05 (dashed line). Experimental points 
taken from Ref. [2]. 



The dispersion curves calculated on the basis of 2 x 2 x 2 
supercell are shown in Fig. [21 They correspond to temper- 
ature T=0. At the same Figure [5] we compare the cal- 
culated dispersion curves with the experimental phonon 
frequencies measured by the inelastic neutron scattering 
[2] . Different symbols of experimental points correspond to 
different phonon polarizations established experimentally. 

The experimental points are about Auj = 0.05w lower 
then the calculated values. One of the reasons for such 
situation could be that experimental values are measured 
at room temperature and anharmonic effects might dem- 
inish the phonon frequencies. Another reason may be re- 
lated with a non-stoichiometric concentration of carbon 
in measured samples, which diminish in average phonon 
frequencies as well. The efect that experimental phonon 
points are at lower frequencies then the calculated ones 
appears as well in other ab-initio calculations pilllj. 

The theory of lattice dynamics provides summation 
rules which follow from translational and rotational invari- 
ances of the crystal. For TiC these invariances lead to two 
equations, which can be used to set the values of on-site 
force constant parameters. If these conditions are satisfied 
the acoustic phonon branches at k=0 would point to w = 
0. Our dispersion curves, without imposed translational- 
rotational invariances point at k= to a; = O.SlTHz. 
Since the invariance conditions should be satisfied within 
the ah initio and direct methods, we have not corrected 
our dispersion curves displayed in Fig. [5] according to this 
effect. 

By sampling the dynamical matrix for many wave vec- 
tors, one can calculate the phonon density function (7(0;), 
and the partial phonon density functions g^^ni^) and 
5a:.c(w). The g{uj) describes the number of phonon fre- 
quencies in an interval around oj, while gx,Ti and gx^c 
specify the number of phonon frequencies in the interval 
around lu but modified by the population of the x displace- 
ment component of either Ti or C atom. The density of 
states functions are shown in Fig. |31 They are convention- 
ally normalized to J g{uj)du! — 1 and J gx.ai^)dui = 1/6, 
where a = Ti or C. The curves show that motions within 
acoustic dispersion curves are almost entirely due to Ti 
atoms. The sum of Titanium density functions {g^^Tii^) + 
gy,Ti{'^) +5z,Ti(w)) fits to the total density function 5(1^) 
below 14 THz. Part of 3(0;) above the gap is mainly due 
to C atoms. Thus, TiC form a rather special crystal, in 
which the Ti atoms, as heavier, form a frame for elas- 
tic motion, and the C atoms vibrate within the optical 
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Fig. 3. Phonon density spectrum g{iij) of the TiC crystal and 
partial phonon density spectra gi^a{ijj), where a = C or Ti and 
i = X, y OT z are the displacement directions. Experimental 
points taken from Ref. [2]. 



modes. We add, however, that the magnitude of the force 
constants between Ti - C and between Ti - Ti and C - C 
at the same distances remains of the same order. Thus, 
TiC does not consist of two weakly bounded subsystems. 
In Fig.[3]the experimentally determined phonon density of 
states [2„26^ is also shown. The agreement within acoustic 
region is quite qood. In optic part the experimental reso- 
lution of « 2.5THz and the nonsoichiometry of carbon in 
the sample lead to broadening of the optic band. 

The elastic constants can be evaluated from the slopes 
of acoustic branches of phonon dispersion curves. Here, 
we have used the invariance conditions to set the acoustic 
modes to w = at k= 0. For the fee structure, the set of 
seven equations for three unknown variables Cij has to be 
solved: pvlj^ = cn; pv^ j, = C44; pvjj^ = Cii+4c44 + 2ci2; 
P'"l,T = Cn + C44 - C12; pv'i^ j^ = cii + 2c44 + C12; pv'i^^Th ^ 
cii — C12; pv^ T3 ~ '^44' where the velocities Vi are the 
slopes of appropriate acoustic branches. Here, a, /? and 
7 indexes describe three main directions in the recipro- 
cal space (0,0,1), (1,1,1) and (0,1,1), respectively. The 
density of TiC p = 4849.6 kg/m^ was found using the cal- 
culated lattice constant. The solution of this set of equa- 
tions, by means of the least-squares fit, results in cn = 5.7, 
C12 = 2.5 and C44 = 1.51 MBar. Comparing these values 
with experimental data (Table [2]) one notices considerable 
discrepancies, especially for C12. Closer look to dispersion 
curves. Fig. [51 assures us that the transverse acoustic mode 
in (1, 1,0) direction and the longitudinal in (1, 1, 1) direc- 
tion do not follow exactly the experimental points. We 
have checked, that such discrepancy of slope is sufficient 
to produce the observed deviations of the calculated elas- 
tic constants, specially of C12. The discrepancy in slopes 
is so small that it could be assigned to the missing force 
constants describing mainly the long range Coulomb in- 
teraction. Indeed, the direct method did not handle well 



interactions exceeding the size of the supercell, except for 
the special points X, L, W mentioned earlier. 



Result 


B 


Cll 


C12 


C44 


Exp. TiCo.91 [3] 


(2.422) 


5.145 


1.060 


1.788 


Exp. TiC [27] 


2.4 


5.00 


1.13 


1.75 


Exp. TiCo.95 [2] 


(2.5) 


5.4 


1.1 


1.8 


Exp. TiCo.89 [H 


(2.5) 


5.2 


1.1 


1.8 


Calculations [8] 


2.2 


4.7 


0.97 


1.67 


Ours, acoustic modes 


(3.57) 


5.7 


2.5 


1.51 


Ours, deform, energy 


2.42 


5.39 


0.94 


1.92 


Ours, stress-strain 


(2.46) 


5.21 


1.09 


1.93 



Table 2. Experimental and calculated bulk modulus B and 
elastic constants dj of the TiC crystal, in units of lO^^Nm"^ = 
IMBar. Values in parenthesis are calculated from B — ^(cn -|- 

2C12). 



To rule out that the discrepancy for C12 does not follow 
from ab-initio calculations, we have derived the bulk mod- 
ulus and the elastic constants by two conventional meth- 
ods [5]. The bulk modulus was calculated in the 1x1x1 
supercell from the total energy as a function of the lattice 
constant. This data have been fitted to a third order poly- 
nomial. Hence, the derivative and the bulk modulus 
B were calculated from the relation: B = —V^, where 
V is the volume of TiC crystallographic unit cell. The B 
values are given in Table [D and they agree quite well with 
experimental data. 

A calculation of cn, C12, C44 elastic constants involved 
a similar procedure, in which the crystal was deformed by 
lengthening and sheer deformations. Deformations from 
0.5% to 3% in length and from 1 to 5 degrees in angle have 
been used. The results of small and large deformations 
are consistent. The deformed lattices have space groups 
Fm3m, I4/mmm and I/mmm, for bulk expansion, elonga- 
tion along z, and sheer modes. In all these deformed lat- 
tices all atoms remain in the high-symmetry sites which 
guaranties an extremum of potential energy at those po- 
sitions. Thus, in the case of cubic TiC, one is entitled 
not to carry on the relaxation of internal degrees of free- 
dom during supercell deformations. The calculated elastic 
constants are compared to experimental data in Table [2] 
They agree with calculated values given in Ref. 8 and are 
in good agreement with experiment. 

We have also checked the above results using stresses 
calculated by CASTEP. The elastic constants were found 
as first derivatives of the stress-strain relationship. For 
that, we have used the same runs of electronic calculations 
as for the elastic constant derivation given above using de- 
formation dependence of energy. The results are included 
in Table [5] and they show slightly better agreement with 
the experimental data. 

In summary, we have shown that the ab initio calcula- 
tions of HF forces, together with the direct method lead 
to satisfactory description of phonons in the TiC crys- 
tal. As usual the frequency of the calculated dispersion 
curves are slightly too high. The acoustic modes are not 
sufficiently accurate, thus the elastic constants calculated 
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from them differ from those found through the supercell 
deformations. The acoustic part of density of states fits 
well to neutron experimental data. 
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